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The local statistical and geometric structure of three-dimensional turbulent flow can be described 
by properties of the velocity gradient tensor. A stochastic model is developed for the Lagrangian 
time evolution of this tensor, in which the exact nonlinear self-stretching term accounts for the 
development of well-known non-Gaussian statistics and geometric alignment trends. The non-local 
pressure and viscous effects are accounted for by a closure that models the material deformation 
history of fluid elements. The resulting stochastic system reproduces many statistical and geometric 
trends observed in numerical and experimental 3D turbulent flows, including anomalous relative 
scaling. 
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Fully developed turbulent flows are omnipresent in the 
natural and man-made environment. Development of 
deeper understanding of fundamental properties of tur- 
bulence is needed for progress in a number of important 
fields such as meteorology, combustion, and astrophysics. 
Despite the highly complex nature of inherently three- 
dimensional velocity fluctuations, turbulent flows exhibit 
universal statistical properties. An example is the k~^^^- 
law of Kolmogorov [l| . Another example is the ubiquity 
of intermittency of longitudinal and transverse Eulerian 
velocity increments between two points 0. Moreover, 
probability density functions (PDFs) of velocity incre- 
ments change with the length-scale between the points. 
Starting from an almost Gaussian density at large scale 
L (i.e. the integral length scale), these PDFs undergo a 
continuous deformation in the inertial range to finish in 
a highly skewed and non Gaussian PDF near the viscous 
scale of turbulence 0, 0. The latter is, equivalent ly, 
also true for the velocity gradients. Recently a simple 
two-equation dynamical system was derived |j that re- 
produces the formation of intermittent tails in the PDFs. 

While much attention has been devoted to the statis- 
tics and anomalous scaling of longitudinal and transverse 
velocity increments, there has been growing interest (see 
e.g. p) in the properties of the full velocity gradient 
tensor Aij = djUi. Aij characterizes variations of all 
velocity components, in all directions. Such additional 
information is required (but unavailable) to model pres- 
sure effects in the system of Refs. 4] and thus to allow 
reproducing stationary statistics. Empirically it has also 
become apparent that A displays a number of interest- 
ing and possibly universal geometric features. For ex- 
ample, the vorticity vector (related to the antisymmetric 
part of A) is preferentially aligned |^ with the eigen- 
vector of the intermediate eigenvalue of the strain-rate 
tensor S = (A + A^)/2, where T stands for transpose. 
Moreover, the preferred state of the local deformation 
is axisymmetric expansion, corresponding to two posi- 
tive and one negative eigenvalues of S. These geomet- 
ric trends have been repeatedly observed in experimental 



and numerical experiments 6], both at the viscous scale 
as well as in the inertial range, for a variety of differ- 
ent flows. These trends can be readily understood from 
the nonlinear self-stretching 0, Q that occurs during the 
Lagrangian evolution of A. However, the resulting so- 
called Restricted Euler (RE) dynamics, obtained by ne- 
glecting viscous diffusion and the non-local anisotropic 
effects of pressure, display unphysical finite-time singu- 
larities. These are due to the absence of regularization 
properties of the neglected viscous and pressure gradi- 
ent terms. Prior models that seek to regularize the RE 
dynamics include a stochastic model in which the nonlin- 
ear term is modified to yield, by construction, log-normal 
statistics of the dissipation 0, a linear damping model 
for the viscous term and the tetrad model in 
which the material deformation history is used to model 
the unclosed pressure Hessian term. Material deforma- 
tion is also tracked in the viscous diffusion closure in Ref. 
[l^ . While each of these models add useful features, a 
model that has no singularities and leads to stationary 
statistics, without tuning the nonlinear term explicitly to 
impose log-normal dissipation statistics, is still lacking. 
The aim of this Letter is to introduce such a model and 
to document its properties. 

The Lagrangian evolution of Aij is governed by the 
gradient of the incompressible Navier- Stokes equations: 



dAij 



— Aik.A, 



ik -^kj 



dxjdxn 



(1) 



where u is the kinematic viscosity, p is the pressure di- 
vided by density, and d/dt the Lagrangian (material) 
derivative. An = at all times. The last two terms in Eq. 
(QJ are unclosed. If the pressure Hessian dfjp is assumed 
to be an isotropic tensor, its trace can be expressed in 
terms of an invariant of A which yields, together with ne- 
glect of the viscous term, to the above-mentioned, closed, 
RE system 8j. Yet, it is well-known that it is unphys- 
ical to assume that dfjp is isotropic, given the complex 
anisotropic effects of pressure gradient. 

If instead we focus on changes of local pressure with 



2 




FIG. 1: Sketch of how an initially uncertain (and thus mod- 
eled as isotropic) material element is mapped onto the present 
position X at time t, reflecting recent deformations. 



changes of past fluid particle locations (X) at some early 
time in the Lagrangian history (i.e. focus on the La- 
grangian pressure Hessian Pmn = d'^p/dX^OXn, where 
p is evaluated at present time t but as function of initial 
positions), the assumption of isotropy is better justified. 
This is based on the idea that any causal relationship 
between the initial time and the present has been lost 
due to the stochastic nature of turbulent dispersion. The 
sketch in Fig. meant to describe how an initially un- 
certain (and thus modelled as isotropic) material shape 
is mapped onto the present location with a deformed 
shape that mirrors the recent local deformations due to 
the velocity gradient history. The notation is as follows: 
x(t) denotes the present position of interest, at time t. 
AltQ^t • X X is the Lagrangian path map which 
gives the Eulerian position x at time t of a fluid particle 
initially located at the position X at time to. By virtue of 
incompressibility, this map is invertible and its Jacobian 
(the deformation gradient tensor), Dj^ = dxi/dXj has 
determinant det(D) = 1 at any time [ij. We denote its 
inverse by D~j^ = dXi/dxj. The tensor dj = DikDjk is 
called the Cauchy-Green tensor which has been studied in 
turbulent flows numerically and experimentally [l5Lll^. 

The relationship between the Eulerian and Lagrangian 
pressure Hessian is obtained by applying twice the change 
of variables d/dxj = {dXm/ dxj)d / dXm^ and neglect- 
ing djdXm. /dxn)/ dxj (i.e. neglecting spatial variations 
of D""*" pj). Then, the main closure hypothesis is that 
the Lagrangian pressure Hessian, is isotropic (i.e. 

Pmn = Pkk^mn/^^ whcrc Smn is the Kroucckcr tensor), 
when the time-delay t — to is long enough to justify loss of 
information. The pressure Hessian can then be rewritten 
according to 



d^p dXm dXn d^p 



dxidxj 



dxi dxj dXmdXn 



(2) 



which could be regarded as a reinterpretation of the 
"tetrad model" 

The dynamics of D are determined by d'D{t)/dt = 
A{t)'D{t). Starting at some initial time from Dij{to) = 
Sij^ the general form of D can be written formally us- 
ing the time-ordered exponential function (exp^), i.e. 



To determine P/e/c, we follow Ref. ^ and use the 

"^P/,/c/3, from 



which Pkk can be solved, leading to 



d'^p 



(3) 

OXiOXj ^qq 

A similar approach can be applied to the viscous 
term, expressing the Laplacian of Aij in the Lagrangian 
frame, as in Eq. Q. The resulting Lagrangian Hessian 
of Aij is modeled by a classical linear damping term, 
namely ud'^ Aij /dXpdXq ^ -SpqAij /{3T). The relax- 
ation time-scale T is chosen to be on the order of the 
integral time-scale. This can be justified by recogniz- 
ing that the distance travelled by a viscous eddy during 
a viscous turn-over or decorrelation time, advected by 
the rms turbulence velocity u\ scales like the Taylor mi- 
croscale, A. Assuming therefore that A is the appropriate 
Lagrangian decorrelation length-scale of A^j, it follows 
that v/{dXY r^vj}? l/T. Finally, the model reads 



d^A, 



1 r-i 

T 3 ' 



(4) 



dXfYidXrYi 

and is reminiscent of mapping closures [l^ . 

Replacing the pressure Hessian and the viscous term in 
Eq. by the modeled terms Eqs. ^ and one can 
show numerically that the finite-time divergence induced 
by the quadratic term is regularized, and each component 
of Aij tends to zero at long times. Next, to generate sta- 
tionary statistics a stochastic forcing term can be added. 
The resulting system, however, is not stationary since 
it depends upon the evolving tensors D and C whose 
time evolutions reflect the non-stationary nature of tur- 
bulent dispersion. For example, on average the largest 
(resp. smallest) eigenvalue of C undergoes exponential 
growth (resp. decrease) in time, whereas the interme- 
diate one remains approximatively constant 0, 0, . 
We remark that in the tetrad model this feature is 
exploited to keep track of changing length scale. Our aim 
here is to develop a statistically stationary description of 
the velocity gradient at a fixed scale (e.g. viscous scale). 

The crucial step of the proposed model is to replace 
the actual slow decorrelation along the Lagrangian tra- 
jectory and the total deformation history [to^t] with a 
perfect correlation of Aij during a time scale r (which 
is thought to be of the order of the Kolmogorov time- 
scale \fyj~^^ where e is the dissipation rate). Correla- 
tions for time-delays longer than r are neglected. It fol- 
lows, using the time-ordered exponential property, that 
D(t) = D(t - r)D^(t), where D^(t) ^ e^^^^\ Further- 
more, we neglect the prior deformation history. Accord- 
ingly, we may define a "stationary Cauchy-Green tensor" 



a(t) = D,(t)D^(t) 



(5) 



Poisson equation V p 



2^ _ 



= G 



qq 



When r decreases (i.e. the Reynolds number 1Ze in- 
creases), at fixed A the restitution strength of the pres- 
sure Hessian model decreases (r = corresponds to an 
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FIG. 2: (a) and (b): PDFs of longitudinal and transverse 
components of the velocity gradient tensor A (normalized by 
its root- mean- square), obtained from time-integration of Eq. 
P for r = 0.2,0.1,0.08 and 0.06. (c) PDF of InAijAij, 
for the same values of F. A Gaussian PDF of unit variance 
(dashed line) is also shown. Arrows direction indicates de- 
creasing F (increasing Tie)- 
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FIG. 3: Contour plots of the logarithm of the joint prob- 
ability of the two invariants of the velocity gradient ten- 
sor (see text) non-dimensionalized by the average strain, i.e. 
Q* = Q/{S^JS^J) and = i?/(5'.,5'^,)^/^ for (a) F = 0.2 and 
(b) F = 0.08. Contours are logarithmically spaced, starting at 
1 and separated by factors of 10. The thick line corresponds 
to zero discriminant (Vieillefosse line). 



isotropic pressure Hessian as in the singular RE sys- 
tem). Without loss of generality, henceforth all variables 
will be scaled with the time-scale T, i.e. t/T ^ t and 
AijT Aij. Combining Eqs. (QJ, Q, © and 
a forcing term, and defining the parameter F = r/T 
(~ Tie ), the following stochastic differential equation 
is finally obtained: 



dA 



Tr(A^) ^ 
Tr(CFi) 



A]dt^dW 



The tensorial noise dW represents physical effects that 
have been neglected, such as action of larger-scale, and 
neighboring, eddies. For simplicity, we assume dW is 
Gaussian and white in time. In the assumed units of 
time, we choose dW = GV2di, where G is a ten- 
sorial Gaussian, delta-correlated noise. Its covariance 
matrix should be consistent with an isotropic, homoge- 
neous, and traceless tensorial field, namely (GijGki) = 
'^hk^ji — \^ij^ki — \^ii^jk nil- WYieii dWij = 0, numerical 
tests show that the finite-time divergence is regularized 
for any initial condition. 

The stochastic differential equation Q is solved nu- 
merically using four different values for F: 0.2, 0.1, 
0.08 and 0.06. A second-order weak predicator-corrector 
scheme is used, with time steps dt = 10~^ {dt = 10~^ 
is used for F = 0.06). Integration times of order 10^ T's 
are used. Time-series of each component of A indi- 
cate stationary behavior. In Figs. El^a-b) we show the 
PDFs of longitudinal (An) and transverse {A12) compo- 
nents for various F values (here and below, all statistics 
are improved by averaging over all available longitudi- 
nal and transverse directions, respectively). When F de- 
creases, velocity gradient PDFs develop slightly longer 
tails. Also, the longitudinal components are negatively 
skewed. 



It has been observed in numerical simulations 
that the pseudo-dissipation Aij Aij is close to lognormal 
for any Reynolds number (as obtained in the station- 
ary diffusion process by specific construction of the 
nonlinear term), and one wonders whether lognormal- 
ity arises in the present model. Fig. (Sljc) presents the 
PDF of the logarithm of the pseudo-dissipation for var- 
ious values of the parameter F. The PDF of hi Aij Aij 
from the model is close (but not exactly equal) to Gaus- 
sian. Note that the finiteness of dissipation implies that 



(6) {Al^)/T'^ = e/{lbiy). It follows that 



is fixed 



through jiyje) 

To further characterize the statistics of A, Fig. 01 
presents the joint PDF of two important invariants of 
A, namely Q = -Tr(A^)/2 and R = -Tr(A^)/3, 
non-dimensionalized by {SijSij). The joint PDF in 
the RQ-plane shows the characteristic teardrope shape 
observed in various numerical and experimental stud- 
ies 0, ^3 is consistent with predominance of 
enstrophy-enstrophy production (top-left quadrant) and 
dissipation-dissipation production (bottom-right quad- 
rant). For decreasing F, the joint PDF becomes more 
elongated along the right tail of the Vieillefosse line, con- 
sistent with data at increasing IZe 0, ^3 • 

Next, the statistics of alignment of the vorticity vector 
= ^ijkAkj with S, and the S-eigenvalues a, (3 and 7 
are quantified. In Fig. El^a) the PDF of cos(6>), where 6 is 
the angle between u; and the S-eigenvector corresponding 
to its intermediate eigenvalue, is shown. Clearly there is 
preferential alignment (as in numerical and experimental 
3D flows 0). To quantify the preferred rate of strain 
state, we display in Fig. El^b) the PDF of the parame- 
ter 5* = -3VQapj/{a'^ 4-/32+ 72)3/2^ As in real flows 
^6], the PDF of s* is shifted towards a peak at s* = 1 
(axisymmetric expansion). 

An important feature of small-scale turbulence is seal- 
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FIG. 4: Alignment trends and preferred strain-rate state. In 
(a) is displayed the PDF of the cosine of the angle between 
vorticity and the intermediate eigenvector of the strain-rate 
tensor, showing preferential alignment. In (b), the PDF of 
the strain-rate state s* is presented. 



ing of higher-order moments with IZe [23 , i-e. |^) 



7^; 



Hp) 



Regular K41 scaling corresponds to J-{p) = p/2 



while deviations indicate anomalous scaling. How- 
ever, the simple assumption to take the forcing term W 
Gaussian and delta correlated in time is expected to be 
realistic at most for a limited range of Reynolds num- 
bers. Therefore, we present results in terms of relative 
scaling which utilizes the above relation for p = 2 to 
obtain IZe ~ (^n) (using F{2) = 1 from the condi- 
tion of finite dissipation), and thus ~ (^n)*^*^^^. 
Shown in Fig. Oare p-order moments of An and ^412, as 
functions of the second-order moments, and varying pa- 
rameter F. Deviations from the dashed lines (K41 case 
with slope p/2) are consistent with anomalous scaling. 
Since PDFs of normalized An and A12 change with F or 
T^e, their statistics cannot follow K41 scaling. The solid 
lines in Fig. [3 use the multifractal formalism: J-{p) = 
- mmh{[p{h-l)^l-V{h)]/{h^l)} and V{h) is the clas- 
sical singularity spectrum 2] . The latter is used here with 
a parabolic approximation T>{h) = 1 — {h — Ci)^/(2c2), 
with ci = 1 /3 + 3c2 /2 0, [l^ , and thus a single unknown 
parameter C2 (c2 = /i/9, where /i is the usual intermit- 
tency exponent). The numerical results can thus be used 
to determine C2 from the model by fitting the slopes in 
Fig. [51 The solid lines are for a parameter C2 = 0.025 (or 
fi = 9c2 ^ 0.22) for the longitudinal, and C2 = 0.040 for 
the transverse cases. These values are in excellent agree- 
ment with values found from experiments and DNS Q^Q- 
The longitudinal derivative skewness factor S shows char- 
acteristic values near —0.5. 

In conclusion, building on several prior works 0,0,0, 
H^ . a new model has been proposed for the anisotropic 
part of the pressure Hessian and the viscous diffusion 
term entering in the Lagrangian evolution equation for 
the velocity gradient tensor A. The system predicts a va- 
riety of local, statistical, geometric and anomalous scal- 
ing properties of 3-D turbulence. Results are obtained 
within a limited range of the parameter F, or Reynolds 





FIG. 5: Relative scaling of velocity gradient moments (o) in 
the (a) longitudinal direction and (b) transverse direction for 
various orders p. Different points are for various F (from left 
to right F = 0.2, 0.15, 0.1, 0.09, 0.08, 0.07, 0.06, and 0.05). In 
(a) the skewness coefficient of longitudinal components S = 
{All) / (AIi)^^"^ is also shown (• using the right scale). Solid 
lines denote predictions from multifractal scaling, dashed lines 
are Kolmogorov (1941) non- anomalous scaling. 



number T^g- When tests are done with F below 0.05, the 
PDFs of velocity increments, of R and Q, and alignment 
trends become less realistic. This is due possibly to the 
limitations imposed by the assumption of Gaussian forc- 
ing. More work is needed to extend the approach to ar- 
bitrarily high Reynolds numbers, possibly by adding ad- 
ditional degrees of freedom to the model or by modifying 
the type of forcing. Moreover, establishing connections 
with the statistics of Lagrangian structure functions (ve- 
locity increments in time instead of the spatial variations 
described by Aij) requires additional models to describe 
jointly the Lagrangian evolution of velocity and velocity 
gradients. 
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